{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# E_vs_r_scan calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The E_vs_r_scan calculation style calculation creates a plot of the cohesive energy vs interatomic spacing, $r$, for a given atomic system. The system size is uniformly scaled ($b/a$ and $c/a$ ratios held fixed) and the energy is calculated at a number of sizes without relaxing the system. All box sizes corresponding to energy minima are identified.\n",
    "\n",
    "This calculation was created as a quick method for scanning the phase space of a crystal structure with a given potential in order to identify starting guesses for further structure refinement calculations.\n",
    "\n",
    "### Version notes\n",
    "\n",
    "- 2018-07-09: Notebook added.\n",
    "- 2019-07-30: Description updated and small changes due to iprPy version.\n",
    "- 2020-05-22: Version 0.10 update - potentials now loaded from database.\n",
    "- 2020-09-22: Setup and parameter definitions streamlined.\n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "- The minima identified by this calculation do not guarantee that the associated crystal structure will be stable as no relaxation is performed by this calculation. Upon relaxation, the atomic positions and box dimensions may transform the system to a different structure.\n",
    "- It is possible that the calculation may miss an existing minima for a crystal structure if it is outside the range of $r$ values scanned, or has $b/a$, $c/a$ values far from the ideal.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "An initial system (and corresponding unit cell system) is supplied. The $r/a$ ratio is identified from the unit cell. The system is then uniformly scaled to all $r_i$ values in the range to be explored and the energy for each is evaluated using LAMMPS and \"run 0\" command, i.e. no relaxations are performed.\n",
    "\n",
    "In identifying energy minima along the curve, only the explored values are used without interpolation. In this way, the possible energy minima structures are identified for $r_i$ where $E(r_i) < E(r_{i-1})$ and $E(r_i) < E(r_{i+1})$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1. Library imports\n",
    "\n",
    "Import libraries needed by the calculation. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2021-03-08 using iprPy version 0.10.4\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "from pathlib import Path\n",
    "import os\n",
    "import datetime\n",
    "from copy import deepcopy\n",
    "from math import floor\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np\n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import additional libraries for plotting. The external libraries used are:\n",
    "\n",
    "- [bokeh](http://bokeh.pydata.org/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bokeh version = 1.2.0\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div class=\"bk-root\">\n",
       "        <a href=\"https://bokeh.pydata.org\" target=\"_blank\" class=\"bk-logo bk-logo-small bk-logo-notebook\"></a>\n",
       "        <span id=\"1001\">Loading BokehJS ...</span>\n",
       "    </div>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "(function(root) {\n",
       "  function now() {\n",
       "    return new Date();\n",
       "  }\n",
       "\n",
       "  var force = true;\n",
       "\n",
       "  if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n",
       "    root._bokeh_onload_callbacks = [];\n",
       "    root._bokeh_is_loading = undefined;\n",
       "  }\n",
       "\n",
       "  var JS_MIME_TYPE = 'application/javascript';\n",
       "  var HTML_MIME_TYPE = 'text/html';\n",
       "  var EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n",
       "  var CLASS_NAME = 'output_bokeh rendered_html';\n",
       "\n",
       "  /**\n",
       "   * Render data to the DOM node\n",
       "   */\n",
       "  function render(props, node) {\n",
       "    var script = document.createElement(\"script\");\n",
       "    node.appendChild(script);\n",
       "  }\n",
       "\n",
       "  /**\n",
       "   * Handle when an output is cleared or removed\n",
       "   */\n",
       "  function handleClearOutput(event, handle) {\n",
       "    var cell = handle.cell;\n",
       "\n",
       "    var id = cell.output_area._bokeh_element_id;\n",
       "    var server_id = cell.output_area._bokeh_server_id;\n",
       "    // Clean up Bokeh references\n",
       "    if (id != null && id in Bokeh.index) {\n",
       "      Bokeh.index[id].model.document.clear();\n",
       "      delete Bokeh.index[id];\n",
       "    }\n",
       "\n",
       "    if (server_id !== undefined) {\n",
       "      // Clean up Bokeh references\n",
       "      var cmd = \"from bokeh.io.state import curstate; print(curstate().uuid_to_server['\" + server_id + \"'].get_sessions()[0].document.roots[0]._id)\";\n",
       "      cell.notebook.kernel.execute(cmd, {\n",
       "        iopub: {\n",
       "          output: function(msg) {\n",
       "            var id = msg.content.text.trim();\n",
       "            if (id in Bokeh.index) {\n",
       "              Bokeh.index[id].model.document.clear();\n",
       "              delete Bokeh.index[id];\n",
       "            }\n",
       "          }\n",
       "        }\n",
       "      });\n",
       "      // Destroy server and session\n",
       "      var cmd = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n",
       "      cell.notebook.kernel.execute(cmd);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  /**\n",
       "   * Handle when a new output is added\n",
       "   */\n",
       "  function handleAddOutput(event, handle) {\n",
       "    var output_area = handle.output_area;\n",
       "    var output = handle.output;\n",
       "\n",
       "    // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n",
       "    if ((output.output_type != \"display_data\") || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n",
       "      return\n",
       "    }\n",
       "\n",
       "    var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n",
       "\n",
       "    if (output.metadata[EXEC_MIME_TYPE][\"id\"] !== undefined) {\n",
       "      toinsert[toinsert.length - 1].firstChild.textContent = output.data[JS_MIME_TYPE];\n",
       "      // store reference to embed id on output_area\n",
       "      output_area._bokeh_element_id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n",
       "    }\n",
       "    if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n",
       "      var bk_div = document.createElement(\"div\");\n",
       "      bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n",
       "      var script_attrs = bk_div.children[0].attributes;\n",
       "      for (var i = 0; i < script_attrs.length; i++) {\n",
       "        toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\n",
       "      }\n",
       "      // store reference to server id on output_area\n",
       "      output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function register_renderer(events, OutputArea) {\n",
       "\n",
       "    function append_mime(data, metadata, element) {\n",
       "      // create a DOM node to render to\n",
       "      var toinsert = this.create_output_subarea(\n",
       "        metadata,\n",
       "        CLASS_NAME,\n",
       "        EXEC_MIME_TYPE\n",
       "      );\n",
       "      this.keyboard_manager.register_events(toinsert);\n",
       "      // Render to node\n",
       "      var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n",
       "      render(props, toinsert[toinsert.length - 1]);\n",
       "      element.append(toinsert);\n",
       "      return toinsert\n",
       "    }\n",
       "\n",
       "    /* Handle when an output is cleared or removed */\n",
       "    events.on('clear_output.CodeCell', handleClearOutput);\n",
       "    events.on('delete.Cell', handleClearOutput);\n",
       "\n",
       "    /* Handle when a new output is added */\n",
       "    events.on('output_added.OutputArea', handleAddOutput);\n",
       "\n",
       "    /**\n",
       "     * Register the mime type and append_mime function with output_area\n",
       "     */\n",
       "    OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n",
       "      /* Is output safe? */\n",
       "      safe: true,\n",
       "      /* Index of renderer in `output_area.display_order` */\n",
       "      index: 0\n",
       "    });\n",
       "  }\n",
       "\n",
       "  // register the mime type if in Jupyter Notebook environment and previously unregistered\n",
       "  if (root.Jupyter !== undefined) {\n",
       "    var events = require('base/js/events');\n",
       "    var OutputArea = require('notebook/js/outputarea').OutputArea;\n",
       "\n",
       "    if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n",
       "      register_renderer(events, OutputArea);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  \n",
       "  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n",
       "    root._bokeh_timeout = Date.now() + 5000;\n",
       "    root._bokeh_failed_load = false;\n",
       "  }\n",
       "\n",
       "  var NB_LOAD_WARNING = {'data': {'text/html':\n",
       "     \"<div style='background-color: #fdd'>\\n\"+\n",
       "     \"<p>\\n\"+\n",
       "     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n",
       "     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n",
       "     \"</p>\\n\"+\n",
       "     \"<ul>\\n\"+\n",
       "     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n",
       "     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n",
       "     \"</ul>\\n\"+\n",
       "     \"<code>\\n\"+\n",
       "     \"from bokeh.resources import INLINE\\n\"+\n",
       "     \"output_notebook(resources=INLINE)\\n\"+\n",
       "     \"</code>\\n\"+\n",
       "     \"</div>\"}};\n",
       "\n",
       "  function display_loaded() {\n",
       "    var el = document.getElementById(\"1001\");\n",
       "    if (el != null) {\n",
       "      el.textContent = \"BokehJS is loading...\";\n",
       "    }\n",
       "    if (root.Bokeh !== undefined) {\n",
       "      if (el != null) {\n",
       "        el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n",
       "      }\n",
       "    } else if (Date.now() < root._bokeh_timeout) {\n",
       "      setTimeout(display_loaded, 100)\n",
       "    }\n",
       "  }\n",
       "\n",
       "\n",
       "  function run_callbacks() {\n",
       "    try {\n",
       "      root._bokeh_onload_callbacks.forEach(function(callback) {\n",
       "        if (callback != null)\n",
       "          callback();\n",
       "      });\n",
       "    } finally {\n",
       "      delete root._bokeh_onload_callbacks\n",
       "    }\n",
       "    console.debug(\"Bokeh: all callbacks have finished\");\n",
       "  }\n",
       "\n",
       "  function load_libs(css_urls, js_urls, callback) {\n",
       "    if (css_urls == null) css_urls = [];\n",
       "    if (js_urls == null) js_urls = [];\n",
       "\n",
       "    root._bokeh_onload_callbacks.push(callback);\n",
       "    if (root._bokeh_is_loading > 0) {\n",
       "      console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n",
       "      return null;\n",
       "    }\n",
       "    if (js_urls == null || js_urls.length === 0) {\n",
       "      run_callbacks();\n",
       "      return null;\n",
       "    }\n",
       "    console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n",
       "    root._bokeh_is_loading = css_urls.length + js_urls.length;\n",
       "\n",
       "    function on_load() {\n",
       "      root._bokeh_is_loading--;\n",
       "      if (root._bokeh_is_loading === 0) {\n",
       "        console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n",
       "        run_callbacks()\n",
       "      }\n",
       "    }\n",
       "\n",
       "    function on_error() {\n",
       "      console.error(\"failed to load \" + url);\n",
       "    }\n",
       "\n",
       "    for (var i = 0; i < css_urls.length; i++) {\n",
       "      var url = css_urls[i];\n",
       "      const element = document.createElement(\"link\");\n",
       "      element.onload = on_load;\n",
       "      element.onerror = on_error;\n",
       "      element.rel = \"stylesheet\";\n",
       "      element.type = \"text/css\";\n",
       "      element.href = url;\n",
       "      console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n",
       "      document.body.appendChild(element);\n",
       "    }\n",
       "\n",
       "    for (var i = 0; i < js_urls.length; i++) {\n",
       "      var url = js_urls[i];\n",
       "      var element = document.createElement('script');\n",
       "      element.onload = on_load;\n",
       "      element.onerror = on_error;\n",
       "      element.async = false;\n",
       "      element.src = url;\n",
       "      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n",
       "      document.head.appendChild(element);\n",
       "    }\n",
       "  };var element = document.getElementById(\"1001\");\n",
       "  if (element == null) {\n",
       "    console.error(\"Bokeh: ERROR: autoload.js configured with elementid '1001' but no matching script tag was found. \")\n",
       "    return false;\n",
       "  }\n",
       "\n",
       "  function inject_raw_css(css) {\n",
       "    const element = document.createElement(\"style\");\n",
       "    element.appendChild(document.createTextNode(css));\n",
       "    document.body.appendChild(element);\n",
       "  }\n",
       "\n",
       "  var js_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-1.2.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-1.2.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-1.2.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-gl-1.2.0.min.js\"];\n",
       "  var css_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-1.2.0.min.css\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-1.2.0.min.css\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-1.2.0.min.css\"];\n",
       "\n",
       "  var inline_js = [\n",
       "    function(Bokeh) {\n",
       "      Bokeh.set_log_level(\"info\");\n",
       "    },\n",
       "    \n",
       "    function(Bokeh) {\n",
       "      \n",
       "    },\n",
       "    function(Bokeh) {} // ensure no trailing comma for IE\n",
       "  ];\n",
       "\n",
       "  function run_inline_js() {\n",
       "    \n",
       "    if ((root.Bokeh !== undefined) || (force === true)) {\n",
       "      for (var i = 0; i < inline_js.length; i++) {\n",
       "        inline_js[i].call(root, root.Bokeh);\n",
       "      }if (force === true) {\n",
       "        display_loaded();\n",
       "      }} else if (Date.now() < root._bokeh_timeout) {\n",
       "      setTimeout(run_inline_js, 100);\n",
       "    } else if (!root._bokeh_failed_load) {\n",
       "      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n",
       "      root._bokeh_failed_load = true;\n",
       "    } else if (force !== true) {\n",
       "      var cell = $(document.getElementById(\"1001\")).parents('.cell').data().cell;\n",
       "      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n",
       "    }\n",
       "\n",
       "  }\n",
       "\n",
       "  if (root._bokeh_is_loading === 0) {\n",
       "    console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n",
       "    run_inline_js();\n",
       "  } else {\n",
       "    load_libs(css_urls, js_urls, function() {\n",
       "      console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n",
       "      run_inline_js();\n",
       "    });\n",
       "  }\n",
       "}(window));"
      ],
      "application/vnd.bokehjs_load.v0+json": "\n(function(root) {\n  function now() {\n    return new Date();\n  }\n\n  var force = true;\n\n  if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n    root._bokeh_onload_callbacks = [];\n    root._bokeh_is_loading = undefined;\n  }\n\n  \n\n  \n  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n    root._bokeh_timeout = Date.now() + 5000;\n    root._bokeh_failed_load = false;\n  }\n\n  var NB_LOAD_WARNING = {'data': {'text/html':\n     \"<div style='background-color: #fdd'>\\n\"+\n     \"<p>\\n\"+\n     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n     \"</p>\\n\"+\n     \"<ul>\\n\"+\n     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n     \"</ul>\\n\"+\n     \"<code>\\n\"+\n     \"from bokeh.resources import INLINE\\n\"+\n     \"output_notebook(resources=INLINE)\\n\"+\n     \"</code>\\n\"+\n     \"</div>\"}};\n\n  function display_loaded() {\n    var el = document.getElementById(\"1001\");\n    if (el != null) {\n      el.textContent = \"BokehJS is loading...\";\n    }\n    if (root.Bokeh !== undefined) {\n      if (el != null) {\n        el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n      }\n    } else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(display_loaded, 100)\n    }\n  }\n\n\n  function run_callbacks() {\n    try {\n      root._bokeh_onload_callbacks.forEach(function(callback) {\n        if (callback != null)\n          callback();\n      });\n    } finally {\n      delete root._bokeh_onload_callbacks\n    }\n    console.debug(\"Bokeh: all callbacks have finished\");\n  }\n\n  function load_libs(css_urls, js_urls, callback) {\n    if (css_urls == null) css_urls = [];\n    if (js_urls == null) js_urls = [];\n\n    root._bokeh_onload_callbacks.push(callback);\n    if (root._bokeh_is_loading > 0) {\n      console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n      return null;\n    }\n    if (js_urls == null || js_urls.length === 0) {\n      run_callbacks();\n      return null;\n    }\n    console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n    root._bokeh_is_loading = css_urls.length + js_urls.length;\n\n    function on_load() {\n      root._bokeh_is_loading--;\n      if (root._bokeh_is_loading === 0) {\n        console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n        run_callbacks()\n      }\n    }\n\n    function on_error() {\n      console.error(\"failed to load \" + url);\n    }\n\n    for (var i = 0; i < css_urls.length; i++) {\n      var url = css_urls[i];\n      const element = document.createElement(\"link\");\n      element.onload = on_load;\n      element.onerror = on_error;\n      element.rel = \"stylesheet\";\n      element.type = \"text/css\";\n      element.href = url;\n      console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n      document.body.appendChild(element);\n    }\n\n    for (var i = 0; i < js_urls.length; i++) {\n      var url = js_urls[i];\n      var element = document.createElement('script');\n      element.onload = on_load;\n      element.onerror = on_error;\n      element.async = false;\n      element.src = url;\n      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n      document.head.appendChild(element);\n    }\n  };var element = document.getElementById(\"1001\");\n  if (element == null) {\n    console.error(\"Bokeh: ERROR: autoload.js configured with elementid '1001' but no matching script tag was found. \")\n    return false;\n  }\n\n  function inject_raw_css(css) {\n    const element = document.createElement(\"style\");\n    element.appendChild(document.createTextNode(css));\n    document.body.appendChild(element);\n  }\n\n  var js_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-1.2.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-1.2.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-1.2.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-gl-1.2.0.min.js\"];\n  var css_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-1.2.0.min.css\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-1.2.0.min.css\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-1.2.0.min.css\"];\n\n  var inline_js = [\n    function(Bokeh) {\n      Bokeh.set_log_level(\"info\");\n    },\n    \n    function(Bokeh) {\n      \n    },\n    function(Bokeh) {} // ensure no trailing comma for IE\n  ];\n\n  function run_inline_js() {\n    \n    if ((root.Bokeh !== undefined) || (force === true)) {\n      for (var i = 0; i < inline_js.length; i++) {\n        inline_js[i].call(root, root.Bokeh);\n      }if (force === true) {\n        display_loaded();\n      }} else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(run_inline_js, 100);\n    } else if (!root._bokeh_failed_load) {\n      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n      root._bokeh_failed_load = true;\n    } else if (force !== true) {\n      var cell = $(document.getElementById(\"1001\")).parents('.cell').data().cell;\n      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n    }\n\n  }\n\n  if (root._bokeh_is_loading === 0) {\n    console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n    run_inline_js();\n  } else {\n    load_libs(css_urls, js_urls, function() {\n      console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n      run_inline_js();\n    });\n  }\n}(window));"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import bokeh\n",
    "print('Bokeh version =', bokeh.__version__)\n",
    "from bokeh.plotting import figure, output_file, show\n",
    "from bokeh.embed import components\n",
    "from bokeh.resources import Resources\n",
    "from bokeh.io import output_notebook\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'E_vs_r_scan'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)\n",
    "    \n",
    "# Initialize connection to library\n",
    "library = iprPy.Library(load=['lammps_potentials'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1. Specify system-specific paths\n",
    "\n",
    "- __lammps_command__ is the LAMMPS command to use (required).\n",
    "\n",
    "- __mpi_command__ MPI command for running LAMMPS in parallel. A value of None will run simulations serially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_command = 'lmp_serial'\n",
    "mpi_command = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Load interatomic potential\n",
    "\n",
    "- __potential_name__ gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.  \n",
    "\n",
    "- __potential__ is an atomman.lammps.Potential object (required)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'\n",
    "\n",
    "# Retrieve potential and parameter file(s)\n",
    "potential = library.get_lammps_potential(id=potential_name, getfiles=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3. Load initial unit cell system\n",
    "\n",
    "- __ucell__ is an atomman.System representing a fundamental unit cell of the system (required).  Here, this is generated using the load parameters and symbols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.520,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  3.520,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  3.520]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.760 |   1.760\n",
      "      2 |       1 |   1.760 |   0.000 |   1.760\n",
      "      3 |       1 |   1.760 |   1.760 |   0.000\n"
     ]
    }
   ],
   "source": [
    "# Create ucell by loading prototype record\n",
    "ucell = am.load('prototype', 'A1--Cu--fcc', symbols='Ni', a=3.52)\n",
    "\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4. Modify system\n",
    "\n",
    "- __sizemults__ list of three integers specifying how many times the ucell vectors of $a$, $b$ and $c$ are replicated in creating system.\n",
    "\n",
    "- __system__ is an atomman.System to perform the scan on (required). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# of atoms in system = 108\n"
     ]
    }
   ],
   "source": [
    "sizemults = [3, 3, 3]\n",
    "\n",
    "# Generate system by supersizing ucell\n",
    "system = ucell.supersize(*sizemults)\n",
    "print('# of atoms in system =', system.natoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5. Specify calculation-specific run parameters\n",
    "\n",
    "- __rmin__ is the minimum r spacing to use.\n",
    "- __rmax__ is the minimum r spacing to use.\n",
    "- __rsteps__ is the number of r spacing steps to evaluate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmin = uc.set_in_units(2.0, 'angstrom')\n",
    "rmax = uc.set_in_units(6.0, 'angstrom')\n",
    "rsteps = 200"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1. run0.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('run0.template', 'w') as f:\n",
    "    f.write(\"\"\"#LAMMPS input script that evaluates a system's energy and pressure without relaxing\n",
    "\n",
    "box tilt large\n",
    "\n",
    "<atomman_system_pair_info>\n",
    "\n",
    "variable peatom equal pe/atoms\n",
    "\n",
    "thermo_style custom step lx ly lz pxx pyy pzz pe v_peatom\n",
    "thermo_modify format float %.13e\n",
    "\n",
    "run 0\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2. e_vs_r()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def e_vs_r(lammps_command, system, potential,\n",
    "           mpi_command=None, ucell=None, \n",
    "           rmin=uc.set_in_units(2.0, 'angstrom'), \n",
    "           rmax=uc.set_in_units(6.0, 'angstrom'), rsteps=200):\n",
    "    \"\"\"\n",
    "    Performs a cohesive energy scan over a range of interatomic spaces, r.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    system : atomman.System\n",
    "        The system to perform the calculation on.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    ucell : atomman.System, optional\n",
    "        The fundamental unit cell correspodning to system.  This is used to\n",
    "        convert system dimensions to cell dimensions. If not given, ucell will\n",
    "        be taken as system.\n",
    "    rmin : float, optional\n",
    "        The minimum r spacing to use (default value is 2.0 angstroms).\n",
    "    rmax : float, optional\n",
    "        The maximum r spacing to use (default value is 6.0 angstroms).\n",
    "    rsteps : int, optional\n",
    "        The number of r spacing steps to evaluate (default value is 200).\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'r_values'** (*numpy.array of float*) - All interatomic spacings,\n",
    "          r, explored.\n",
    "        - **'a_values'** (*numpy.array of float*) - All unit cell a lattice\n",
    "          constants corresponding to the values explored.\n",
    "        - **'Ecoh_values'** (*numpy.array of float*) - The computed cohesive\n",
    "          energies for each r value.\n",
    "        - **'min_cell'** (*list of atomman.System*) - Systems corresponding to\n",
    "          the minima identified in the Ecoh_values.\n",
    "    \"\"\"\n",
    "    # Build filedict if function was called from iprPy\n",
    "    try:\n",
    "        assert __name__ == pkg_name\n",
    "        calc = iprPy.load_calculation(calculation_style)\n",
    "        filedict = calc.filedict\n",
    "    except:\n",
    "        filedict = {}\n",
    "\n",
    "    # Make system a deepcopy of itself (protect original from changes)\n",
    "    system = deepcopy(system)\n",
    "    \n",
    "    # Set ucell = system if ucell not given\n",
    "    if ucell is None:\n",
    "        ucell = system\n",
    "    \n",
    "    # Calculate the r/a ratio for the unit cell\n",
    "    r_a = r_a_ratio(ucell)\n",
    "    \n",
    "    # Get ratios of lx, ly, and lz of system relative to a of ucell\n",
    "    lx_a = system.box.a / ucell.box.a\n",
    "    ly_a = system.box.b / ucell.box.a\n",
    "    lz_a = system.box.c / ucell.box.a\n",
    "    alpha = system.box.alpha\n",
    "    beta =  system.box.beta\n",
    "    gamma = system.box.gamma\n",
    " \n",
    "    # Build lists of values\n",
    "    r_values = np.linspace(rmin, rmax, rsteps)\n",
    "    a_values = r_values / r_a\n",
    "    Ecoh_values = np.empty(rsteps)\n",
    "    \n",
    "    # Loop over values\n",
    "    for i in range(rsteps):\n",
    "        \n",
    "        # Rescale system's box\n",
    "        a = a_values[i]\n",
    "        system.box_set(a = a * lx_a, \n",
    "                       b = a * ly_a, \n",
    "                       c = a * lz_a, \n",
    "                       alpha=alpha, beta=beta, gamma=gamma, scale=True)\n",
    "        \n",
    "        # Get lammps units\n",
    "        lammps_units = lmp.style.unit(potential.units)\n",
    "        \n",
    "        # Define lammps variables\n",
    "        lammps_variables = {}\n",
    "        system_info = system.dump('atom_data', f='atom.dat',\n",
    "                                  potential=potential,\n",
    "                                  return_pair_info=True)\n",
    "        lammps_variables['atomman_system_pair_info'] = system_info\n",
    "        \n",
    "        # Write lammps input script\n",
    "        template_file = 'run0.template'\n",
    "        lammps_script = 'run0.in'\n",
    "        template = iprPy.tools.read_calc_file(template_file, filedict)\n",
    "        with open(lammps_script, 'w') as f:\n",
    "            f.write(iprPy.tools.filltemplate(template, lammps_variables,\n",
    "                                             '<', '>'))\n",
    "        \n",
    "        # Run lammps and extract data\n",
    "        try:\n",
    "            output = lmp.run(lammps_command, lammps_script, mpi_command)\n",
    "        except:\n",
    "            Ecoh_values[i] = np.nan\n",
    "        else:\n",
    "            thermo = output.simulations[0]['thermo']\n",
    "            \n",
    "            if output.lammps_date < datetime.date(2016, 8, 1):\n",
    "                Ecoh_values[i] = uc.set_in_units(thermo.peatom.values[-1],\n",
    "                                                lammps_units['energy'])\n",
    "            else:\n",
    "                Ecoh_values[i] = uc.set_in_units(thermo.v_peatom.values[-1],\n",
    "                                                lammps_units['energy'])\n",
    "        \n",
    "        # Rename log.lammps\n",
    "        try:\n",
    "            shutil.move('log.lammps', 'run0-'+str(i)+'-log.lammps')\n",
    "        except:\n",
    "            pass\n",
    "\n",
    "    if len(Ecoh_values[np.isfinite(Ecoh_values)]) == 0:\n",
    "        raise ValueError('All LAMMPS runs failed. Potential likely invalid or incompatible.')  \n",
    "    \n",
    "    # Find unit cell systems at the energy minimums\n",
    "    min_cells = []\n",
    "    for i in range(1, rsteps - 1):\n",
    "        if (Ecoh_values[i] < Ecoh_values[i-1]\n",
    "            and Ecoh_values[i] < Ecoh_values[i+1]):\n",
    "            a = a_values[i]\n",
    "            cell = deepcopy(ucell)\n",
    "            cell.box_set(a = a,\n",
    "                         b = a * ucell.box.b / ucell.box.a,\n",
    "                         c = a * ucell.box.c / ucell.box.a, \n",
    "                         alpha=alpha, beta=beta, gamma=gamma, scale=True)\n",
    "            min_cells.append(cell)\n",
    "    \n",
    "    # Collect results\n",
    "    results_dict = {}\n",
    "    results_dict['r_values'] = r_values\n",
    "    results_dict['a_values'] = a_values\n",
    "    results_dict['Ecoh_values'] = Ecoh_values\n",
    "    results_dict['min_cell'] = min_cells\n",
    "    \n",
    "    return results_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.3. r_a_ratio()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def r_a_ratio(ucell):\n",
    "    \"\"\"\n",
    "    Calculates the r/a ratio by identifying the shortest interatomic spacing, r,\n",
    "    for a unit cell.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    ucell : atomman.System\n",
    "        The unit cell system to evaluate.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        The shortest interatomic spacing, r, divided by the unit cell's a\n",
    "        lattice parameter.\n",
    "    \"\"\"\n",
    "    r_a = ucell.box.a\n",
    "    for i in range(ucell.natoms):\n",
    "        for j in range(i):\n",
    "            dmag = np.linalg.norm(ucell.dvect(i,j))\n",
    "            if dmag < r_a:\n",
    "                r_a = dmag\n",
    "    return r_a / ucell.box.a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = e_vs_r(lammps_command, system, potential,\n",
    "                      mpi_command = mpi_command, \n",
    "                      ucell = ucell, \n",
    "                      rmin = rmin, \n",
    "                      rmax = rmax, \n",
    "                      rsteps = rsteps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['r_values', 'a_values', 'Ecoh_values', 'min_cell'])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1. Define units for outputting values\n",
    "\n",
    "- __length_unit__ is the unit of length to display values in.\n",
    "- __energy_unit__ is the unit of energy to display values in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "length_unit = 'angstrom'\n",
    "energy_unit = 'eV'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2. Plot Ecoh vs r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "  <div class=\"bk-root\" id=\"5d76c6d2-14e0-4586-a3a3-3eb92287e555\" data-root-id=\"1002\"></div>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "(function(root) {\n",
       "  function embed_document(root) {\n",
       "    \n",
       "  var docs_json = {\"ee84efb3-958d-4d1e-8022-63559d43405a\":{\"roots\":{\"references\":[{\"attributes\":{\"below\":[{\"id\":\"1013\",\"type\":\"LinearAxis\"}],\"center\":[{\"id\":\"1017\",\"type\":\"Grid\"},{\"id\":\"1022\",\"type\":\"Grid\"}],\"left\":[{\"id\":\"1018\",\"type\":\"LinearAxis\"}],\"plot_width\":800,\"renderers\":[{\"id\":\"1039\",\"type\":\"GlyphRenderer\"}],\"title\":{\"id\":\"1003\",\"type\":\"Title\"},\"toolbar\":{\"id\":\"1029\",\"type\":\"Toolbar\"},\"x_range\":{\"id\":\"1005\",\"type\":\"Range1d\"},\"x_scale\":{\"id\":\"1009\",\"type\":\"LinearScale\"},\"y_range\":{\"id\":\"1007\",\"type\":\"Range1d\"},\"y_scale\":{\"id\":\"1011\",\"type\":\"LinearScale\"}},\"id\":\"1002\",\"subtype\":\"Figure\",\"type\":\"Plot\"},{\"attributes\":{\"callback\":null,\"end\":6.0,\"start\":2.0},\"id\":\"1005\",\"type\":\"Range1d\"},{\"attributes\":{\"callback\":null,\"end\":0,\"start\":-5},\"id\":\"1007\",\"type\":\"Range1d\"},{\"attributes\":{\"ticker\":{\"id\":\"1014\",\"type\":\"BasicTicker\"}},\"id\":\"1017\",\"type\":\"Grid\"},{\"attributes\":{\"callback\":null,\"data\":{\"x\":{\"__ndarray__\":\"AAAAAAAAAEChZcdzKikAQELLjudUUgBA4jBWW397AECDlh3PqaQAQCT85ELUzQBAxWGstv72AEBmx3MqKSABQAYtO55TSQFAp5ICEn5yAUBI+MmFqJsBQOldkfnSxAFAisNYbf3tAUAqKSDhJxcCQMuO51RSQAJAbPSuyHxpAkANWnY8p5ICQK6/PbDRuwJATiUFJPzkAkDvisyXJg4DQJDwkwtRNwNAMVZbf3tgA0DSuyLzpYkDQHIh6mbQsgNAE4ex2vrbA0C07HhOJQUEQFVSQMJPLgRA9rcHNnpXBECWHc+ppIAEQDeDlh3PqQRA2OhdkfnSBEB5TiUFJPwEQBq07HhOJQVAuhm07HhOBUBbf3tgo3cFQPzkQtTNoAVAnUoKSPjJBUA+sNG7IvMFQN4VmS9NHAZAgHtgo3dFBkAg4ScXom4GQMFG74rMlwZAYqy2/vbABkACEn5yIeoGQKR3ReZLEwdARN0MWnY8B0DlQtTNoGUHQIaom0HLjgdAJg5jtfW3B0DIcyopIOEHQGjZ8ZxKCghACT+5EHUzCECqpICEn1wIQEoKSPjJhQhA7G8PbPSuCECM1dbfHtgIQC07nlNJAQlAzqBlx3MqCUBuBi07nlMJQBBs9K7IfAlAsNG7IvOlCUBRN4OWHc8JQPKcSgpI+AlAkgISfnIhCkA0aNnxnEoKQNTNoGXHcwpAdTNo2fGcCkAWmS9NHMYKQLb+9sBG7wpAWGS+NHEYC0D4yYWom0ELQJkvTRzGagtAOpUUkPCTC0Da+tsDG70LQHxgo3dF5gtAHMZq628PDEC9KzJfmjgMQF6R+dLEYQxA//bARu+KDECgXIi6GbQMQEDCTy5E3QxA4ScXom4GDUCCjd4VmS8NQCPzpYnDWA1AxFht/e2BDUBkvjRxGKsNQAUk/ORC1A1AponDWG39DUBH74rMlyYOQOhUUkDCTw5AiLoZtOx4DkApIOEnF6IOQMqFqJtByw5Aa+tvD2z0DkAMUTeDlh0PQKy2/vbARg9ATRzGautvD0DugY3eFZkPQI/nVFJAwg9AME0cxmrrD0Bo2fGcSgoQQDiM1dbfHhBACT+5EHUzEEDa8ZxKCkgQQKqkgISfXBBAeldkvjRxEEBKCkj4yYUQQBu9KzJfmhBA7G8PbPSuEEC8IvOlicMQQIzV1t8e2BBAXIi6GbTsEEAtO55TSQERQP7tgY3eFRFAzqBlx3MqEUCeU0kBCT8RQG4GLTueUxFAP7kQdTNoEUAQbPSuyHwRQOAe2OhdkRFAsNG7IvOlEUCAhJ9ciLoRQFE3g5YdzxFAIupm0LLjEUDynEoKSPgRQMJPLkTdDBJAkgISfnIhEkBjtfW3BzYSQDRo2fGcShJABBu9KzJfEkDUzaBlx3MSQKSAhJ9ciBJAdTNo2fGcEkBG5ksTh7ESQBaZL00cxhJA5ksTh7HaEkC2/vbARu8SQIex2vrbAxNAWGS+NHEYE0AoF6JuBi0TQPjJhaibQRNAyHxp4jBWE0CZL00cxmoTQGriMFZbfxNAOpUUkPCTE0AKSPjJhagTQNr62wMbvRNAq62/PbDRE0B8YKN3ReYTQEwTh7Ha+hNAHMZq628PFEDteE4lBSQUQL0rMl+aOBRAjt4VmS9NFEBekfnSxGEUQC5E3QxadhRA//bARu+KFEDPqaSAhJ8UQKBciLoZtBRAcA9s9K7IFEBAwk8uRN0UQBF1M2jZ8RRA4ScXom4GFUCy2vrbAxsVQIKN3hWZLxVAUkDCTy5EFUAj86WJw1gVQPOlicNYbRVAxFht/e2BFUCUC1E3g5YVQGS+NHEYqxVANXEYq62/FUAFJPzkQtQVQNbW3x7Y6BVAponDWG39FUB2PKeSAhIWQEfvisyXJhZAF6JuBi07FkDoVFJAwk8WQLgHNnpXZBZAiLoZtOx4FkBZbf3tgY0WQCkg4ScXohZA+tLEYay2FkDKhaibQcsWQJo4jNXW3xZAa+tvD2z0FkA7nlNJAQkXQAxRN4OWHRdA3AMbvSsyF0Cstv72wEYXQH1p4jBWWxdATRzGautvF0Aez6mkgIQXQO6Bjd4VmRdAvjRxGKutF0CP51RSQMIXQGCaOIzV1hdAME0cxmrrF0AAAAAAAAAYQA==\",\"dtype\":\"float64\",\"shape\":[200]},\"y\":{\"__ndarray__\":\"NDsO73P07L8z1EPEFifxv9vQaVy1fvS/ETK/ZJ9L+L8kKt3LvVz8v7CNJilvPgDAt48u5/9AAsDR08IWOi8EwCxeCKSpCAbAL6+UKSjMB8AEYbzlDHgJwHC3zpTKCQvA/Yeob2B+DMASpdlXktENwNI5t6jb/w7AwyhTtJIDEMDVX6PkV3MQwFdDcbyuzxDAJUSSDuIZEcARWADZA1QRwFHQf0tbgBHAJfNLoPCgEcDHHFZ3sbcRwPWWR4oMxhHA9yqiymXMEcBSkAEwGcsRwJKi8RQBwhHA2hWq/dOwEcClz2LfyZcRwM1f+X+cdxHAFuv7sWZREcDY0h/0KSYRwN31zV6Y9hDAJQvbrITDEMAFttZs4Y0QwKkfiSCmVhDALRIEbKUeEMDmlvmWvcwPwGQGgPaUWw/AKrDDZrfoDsDuQaHWMnMOwFpIUT0j+g3ArvuumdV8DcBW8rja8PoMwC4qcVZydAzAnglARpnpC8CVC4WUqVoLwBHlkVWdxwrAc5NIIX0wCsB18o5gl5UJwM3lQnBJ9wjAb/t5tzpUCMDK3vMX2KoHwHGBqblN+gbA0rludRNCBsBkbv1I6oEFwBc7i9IQugTAVMplIvnqA8AuNpp5AxUDwFgOuwJPOALAPK2rpIFVAcC8jwoXv28AwMJ94vLBE/+/zPu/FWBM/b8ciCh7To37v3qTUpUr2vm/J2Y+j1M2+L/7P3Yj2qT2v5KdBCVVKPW/MyqfVNPB87+E5lc/YXHyv9zCute0NvG/iOul0kIR8L/M+xKmmwDuv0TlYGbmBey/ts6JSmsw6r/RfUQlDH7ov6qQx++V7Oa/2G4qE8155b8LTOXfdyPkvz9mK7tn5+K/S2OUeIHD4b8rq967w7Xgv2IKx7GdeN+/MiOj8N2q3b8VGUIxTf/bv2fQiulnc9q/sGhVBSAF2b8pGABg1LLXv4ZdoBE3e9a/nrzaeCBd1b/scv2gUlfUv58ErVs4aNO/zOCIG/SN0r8xVw0ci8bRv57pDyhnD9G/ghZlGxFm0L+5jWO7BpHPvywZELlmac6/gMHdBXZRzb9E/+hh80XMvwdD8kVERMu/B9QmMXJKyr+IBGDgE1fJv76RDXYraci/kBfgJf1/x7+b8gbWSprGvyrjkPfitcW/XfT1e3/QxL9GbOpx5ufDv4MNFLot+8K/tYKXN+4Kwr+eSMNI1xfBvzF3KlGZIsC/8rGuistXvr+Bs/JF3mi8v/DHCuDRebq/NtvxZA+MuL8cKV/ZAaG2v+U8YCcWurS/mrNTBrvYsr/1VTbiYP6wv4hFWo/zWK6/YpeDcfLIqr8Pqc8TqE+nv0mznYAA8KO/8kX+JOmsoL+wiIUgkxKbv5w/NKAMEJW/+CAGWRiwjr/kd/SJxOCDv3/AvA+Ee3O/Kft2N0itJL85WHoegqlwPwKIb/oFLYA/yHLRIIkqhz+BOwamx0GNP2MEcPqrM5E/bWXbNN1Lkz9pY76k0PCUP7P3/3wzK5Y/hwVuHMIDlz9kQE04SYOXPxmzxo+mspc/Ue+i+8ialz/hYASor0SXP826VyRouZY/LCGm6woClj84Rl3ZtSeVP1ai2QuEM5Q/TiMfYYIukz+iYD7BnyGSPz1zUzKYFZE/Yy9vmdkSkD+wS4UUwkKOPwpfQw8fkYw/xsZvzeodiz+6ZkzuJ/SJP4o8Qx0sHIk/gkaY6wKbiD+yWYr5x3GIP9LijQ7znIg/Belukd8TiT9dvNunk8iJP8ClHhfEp4o/qCxrXE2Ziz/dWgs3P4GMP/DDyleIQY0/3beQKjy8jT8xR20+T9aNP41opjaHeo0/EuuCxEGcjD/HMtkbrDmLP45RVRAQXYk/d82ARP4chz9rELcRTZuEP/7qw4sY/IE/WKfDQemnfj+ePNFOhG15P4duU8xrcnQ/Etsz9qm5bz/UzayFHphnPy0+xVI2rGA/U+jKV9wZVj9eKCbJh+9KP5JpzqX1Mz0/hVqIX6aEKj8tpja8TbsRP/0zlpsNLOg+ykD5y8+6dT7C///rpAbBvcL//+ukBsG9wv//66QGwb3C///rpAbBvcL//+ukBsG9wv//66QGwb3C///rpAbBvcL//+ukBsG9wv//66QGwb3C///rpAbBvQ==\",\"dtype\":\"float64\",\"shape\":[200]}},\"selected\":{\"id\":\"1047\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"1048\",\"type\":\"UnionRenderers\"}},\"id\":\"1036\",\"type\":\"ColumnDataSource\"},{\"attributes\":{\"line_color\":\"#1f77b4\",\"line_width\":2,\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"1037\",\"type\":\"Line\"},{\"attributes\":{\"overlay\":{\"id\":\"1046\",\"type\":\"BoxAnnotation\"}},\"id\":\"1025\",\"type\":\"BoxZoomTool\"},{\"attributes\":{},\"id\":\"1023\",\"type\":\"PanTool\"},{\"attributes\":{},\"id\":\"1027\",\"type\":\"ResetTool\"},{\"attributes\":{},\"id\":\"1044\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{\"bottom_units\":\"screen\",\"fill_alpha\":{\"value\":0.5},\"fill_color\":{\"value\":\"lightgrey\"},\"left_units\":\"screen\",\"level\":\"overlay\",\"line_alpha\":{\"value\":1.0},\"line_color\":{\"value\":\"black\"},\"line_dash\":[4,4],\"line_width\":{\"value\":2},\"render_mode\":\"css\",\"right_units\":\"screen\",\"top_units\":\"screen\"},\"id\":\"1046\",\"type\":\"BoxAnnotation\"},{\"attributes\":{},\"id\":\"1014\",\"type\":\"BasicTicker\"},{\"attributes\":{\"text\":\"Cohesive Energy vs. Interatomic Spacing\"},\"id\":\"1003\",\"type\":\"Title\"},{\"attributes\":{},\"id\":\"1028\",\"type\":\"HelpTool\"},{\"attributes\":{},\"id\":\"1019\",\"type\":\"BasicTicker\"},{\"attributes\":{\"axis_label\":\"Cohesive Energy (eV/atom)\",\"formatter\":{\"id\":\"1042\",\"type\":\"BasicTickFormatter\"},\"ticker\":{\"id\":\"1019\",\"type\":\"BasicTicker\"}},\"id\":\"1018\",\"type\":\"LinearAxis\"},{\"attributes\":{\"active_drag\":\"auto\",\"active_inspect\":\"auto\",\"active_multi\":null,\"active_scroll\":\"auto\",\"active_tap\":\"auto\",\"tools\":[{\"id\":\"1023\",\"type\":\"PanTool\"},{\"id\":\"1024\",\"type\":\"WheelZoomTool\"},{\"id\":\"1025\",\"type\":\"BoxZoomTool\"},{\"id\":\"1026\",\"type\":\"SaveTool\"},{\"id\":\"1027\",\"type\":\"ResetTool\"},{\"id\":\"1028\",\"type\":\"HelpTool\"}]},\"id\":\"1029\",\"type\":\"Toolbar\"},{\"attributes\":{\"data_source\":{\"id\":\"1036\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"1037\",\"type\":\"Line\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"1038\",\"type\":\"Line\"},\"selection_glyph\":null,\"view\":{\"id\":\"1040\",\"type\":\"CDSView\"}},\"id\":\"1039\",\"type\":\"GlyphRenderer\"},{\"attributes\":{\"source\":{\"id\":\"1036\",\"type\":\"ColumnDataSource\"}},\"id\":\"1040\",\"type\":\"CDSView\"},{\"attributes\":{},\"id\":\"1011\",\"type\":\"LinearScale\"},{\"attributes\":{},\"id\":\"1042\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{},\"id\":\"1026\",\"type\":\"SaveTool\"},{\"attributes\":{\"dimension\":1,\"ticker\":{\"id\":\"1019\",\"type\":\"BasicTicker\"}},\"id\":\"1022\",\"type\":\"Grid\"},{\"attributes\":{},\"id\":\"1048\",\"type\":\"UnionRenderers\"},{\"attributes\":{\"axis_label\":\"r (angstrom)\",\"formatter\":{\"id\":\"1044\",\"type\":\"BasicTickFormatter\"},\"ticker\":{\"id\":\"1014\",\"type\":\"BasicTicker\"}},\"id\":\"1013\",\"type\":\"LinearAxis\"},{\"attributes\":{\"line_alpha\":0.1,\"line_color\":\"#1f77b4\",\"line_width\":2,\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"1038\",\"type\":\"Line\"},{\"attributes\":{},\"id\":\"1009\",\"type\":\"LinearScale\"},{\"attributes\":{},\"id\":\"1024\",\"type\":\"WheelZoomTool\"},{\"attributes\":{},\"id\":\"1047\",\"type\":\"Selection\"}],\"root_ids\":[\"1002\"]},\"title\":\"Bokeh Application\",\"version\":\"1.2.0\"}};\n",
       "  var render_items = [{\"docid\":\"ee84efb3-958d-4d1e-8022-63559d43405a\",\"roots\":{\"1002\":\"5d76c6d2-14e0-4586-a3a3-3eb92287e555\"}}];\n",
       "  root.Bokeh.embed.embed_items_notebook(docs_json, render_items);\n",
       "\n",
       "  }\n",
       "  if (root.Bokeh !== undefined) {\n",
       "    embed_document(root);\n",
       "  } else {\n",
       "    var attempts = 0;\n",
       "    var timer = setInterval(function(root) {\n",
       "      if (root.Bokeh !== undefined) {\n",
       "        embed_document(root);\n",
       "        clearInterval(timer);\n",
       "      }\n",
       "      attempts++;\n",
       "      if (attempts > 100) {\n",
       "        console.log(\"Bokeh: ERROR: Unable to run BokehJS code because BokehJS library is missing\");\n",
       "        clearInterval(timer);\n",
       "      }\n",
       "    }, 10, root)\n",
       "  }\n",
       "})(window);"
      ],
      "application/vnd.bokehjs_exec.v0+json": ""
     },
     "metadata": {
      "application/vnd.bokehjs_exec.v0+json": {
       "id": "1002"
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Ecoh = uc.get_in_units(results_dict['Ecoh_values'], energy_unit)\n",
    "r = uc.get_in_units(results_dict['r_values'], length_unit)\n",
    "\n",
    "Emin = floor(Ecoh.min())\n",
    "if Emin < -10: \n",
    "    Emin = -10\n",
    "    \n",
    "plot = figure(title = f'Cohesive Energy vs. Interatomic Spacing',\n",
    "              plot_width = 800,\n",
    "              plot_height = 600,\n",
    "              x_range = [uc.get_in_units(rmin, 'angstrom'), uc.get_in_units(rmax, 'angstrom')],\n",
    "              y_range = [Emin, 0],              \n",
    "              x_axis_label=f'r ({length_unit})', \n",
    "              y_axis_label=f'Cohesive Energy ({energy_unit}/atom)')\n",
    "\n",
    "plot.line(r, Ecoh, line_width = 2)              \n",
    "\n",
    "show(plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3. List minimum energy configurations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Possible minimum near:\n",
      "a = 3.510660803076929 angstrom\n",
      "b = 3.510660803076929 angstrom\n",
      "c = 3.510660803076929 angstrom\n",
      "\n",
      "Possible minimum near:\n",
      "a = 7.376651646951118 angstrom\n",
      "b = 7.376651646951118 angstrom\n",
      "c = 7.376651646951118 angstrom\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for mincell in results_dict['min_cell']:\n",
    "    print('Possible minimum near:')\n",
    "    print('a =', uc.get_in_units(mincell.box.a, length_unit), length_unit)\n",
    "    print('b =', uc.get_in_units(mincell.box.b, length_unit), length_unit)\n",
    "    print('c =', uc.get_in_units(mincell.box.c, length_unit), length_unit)\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
